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ABSTRACT 

For advection schemes based on fluctuation splitting, a design criterion of optimising the 
time step leads to linear schemes that coincide with those designed for least truncation error. 
A further stage of optimising the time step using a non-linear positivity criterion, leads to 
considerable further gains in resolution. 
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tract No. NAS1-18605 while the author was in residence at the Institute for Computer Applications in 
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1. Introduction 


This note is a contribution to the development of high resolution advection schemes for 
use on (possibly) unstructured triangular meshes. The search is for a scheme that is positive 
(admitting no new extrema) and as accurate as possible. Such schemes have an inherently 
nonlinear dependence on the data, and revert near discontinuities to an essentially first-order 
behavior. It is of interest to know what is the optimum positive first-order scheme, so that 
performance near discontinuities is not excessively degraded. 

The schemes considered in this paper belong to a class introduced in [1], and explored 
more thoroughly in [2]. These papers, and also [3], should be consulted for more background. 

2. Previous Work 

Consider a triangular grid of which Figure 1 shows a generic element. On such a grid, 
we wish to solve the two-dimensional advection equation 


u t + a • Vu = 0, 


( 1 ) 


where & is a constant vector (a, b) T . 

For each triangular element T define the “fluctuation” 

<fi T = J J u t dxdy. 

By Equation (1) and Gauss’ theorem 

<f> T = £ u a • dn. 

If u is assumed to vary linearly over the triangle 

where 1 

K = 2 ~' 

and rii is the inward (scaled) normal to the side opposite Note that 

x> =° 

so that either one or two k t are positive. Geometrically, this corresponds to the vector a 
crossing either one or two sides in the inward direction. 

In [1,2] it is shown that part of a positive update process is to add, in the case of one inflow 
side (which we call a type I triangle, Figure 2a) a quantity A t<f>x to the single downstream 
vertex. If this is Vi, then 

5i< +1 = Siu* + A tf + (•) (4) 

where (•) signifies the contributions, if any, coming from other triangles neighboring v x , and 
Si is one-third the area of those triangles having v x as a vertex. 
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In the difficult case where there are two inflow sides, say those opposite vertices i and j 
(a type II triangle) the proposed scheme is 

£< +1 = + aiAt<f> T + (•) 

Siuf 1 = S jU “ + aj At<f> T + (•) 
where, from a conservation argument, 



CHi -j- CLj = 1. 


It was proved in [1] (see also [2]) that no choice of a», a,- depending only on the geometry 
of the triangle and its relation to the flow vector can give a positive scheme. The nonlinear 
choice, which we call SI, 

oci = K(ui - u k )l<j> T ) 


( 6 ) 


OLj = k 3 (uj - u k )/<f> T J 

was advocated in [1] on the grounds that it permitted a larger timestep than any other choice 
(within a certain class). This scheme was found to be superior to regular finite- volume 
upwinding in [2, 3] and we will point out in the next section that it bears an interesting 
relationship to a recently published scheme of Sidilkover [4]. 


3. A Property of the Scheme SI 

Substituting (6) into (5) gives the simple update method 


5ju“ +1 = SiU J* - kiAt(ui - u k ) -(- (•) 
Sjuf 1 = SjU* - kjAt(u” - v%) + (•) ; 


( 7 ) 


Let us see what (7) reduces to on a square mesh of side h. We assume that a, b are both 
positive with a > 6, and that the squares are divided into triangles along the diagonal with 
positive slope, as in Figure 3a. Freedom to cut the squares along either diagonal greatly 
enhances the performance of the schemes (R. Struijs, H. Deconinck, private communication, 
von Karman Institute). 

Each point then receives contributions from three triangles, A , 5, (7, in which the fluctu- 
ations are 

4> A = ~^ [(q - b)(uo - tt 3 ) fe + - u 3)^] 

4> B = -^[ a(u 0 - u 3 ) h b(u 3 - u 4 ) h] 

<t> c = ~7j[(a - &)(u 5 - u 4 )h + b(u Q - u 4 ) hj. 

In each case, the fluctuations have been split as they would be by SI, and those terms 
contributing to the update of u 0 have been underlined. Summing those terms, we find 

< +1 = «o - ^rK«o - u 3 ) + b(u 3 - I*)]. (8) 
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This is the scheme recently shown by Sidilkover [4] to have the least truncation error 
among the linear monotone schemes for solving (1). It comes in eight different flavors 
depending on the flow direction (Figure 4). (It may be remarked that whether the scheme 
is linear or nonlinear depends on the viewpoint. Clearly (7) has constant coefficients, but 
they change as the vector a moves into different octants. Written in the form (5), (6) the 
expressions that accomplish this switching give the appearance of nonlinearity.) Sidilkover 
[4] calls (7) the jV-scheme because it has a narrow stencil and, not by coincidence, yields 
narrow discontinuities. 

Now consider the squares divided along the other family of diagonals as in Figure 3b. 
The fluctuations in the triangles adjoining 0 are 

(j> A = —ajup — u 3 ) — b(u 2 — u 3 ) 


<f> B = - ajup - u 3 ) - 6(u 0 - u 5 ) 

4> C = —a(u 6 — u 5 ) — bjup — u 5 ) 

and scheme SI will credit the underlined terms to 0. The update formula is 

u o +1 = u o ~ -^~[ a K) “ 7X3 ) + K u o “ «.)] ( 9 ) 

which is just regular upwinding. For steady flows, the dissipation coefficient of (9) is 

■^1 sin 8 cos 0|{| sin 8\ 4- | cos 0|} 

Cm 

whereas the dissipation coefficient of (8) is 

^|sin$cos<?| | cos 9 — sin B\ 

which is, on average, about four times smaller. There seems no reason why one should not 
take advantage of this by choosing the triangulation appropriate to the problem. In the 
non-linear case, of course, the choice might change with time. 


4. Type II Triangles Revisited 

The observation is made in [1, 2] that within each triangle T the advection velocity a can 
be replaced by 

a T = a + AaJ (10) 

where aj is orthogonal to Vu T and A is an arbitrary parameter. This is tantamount to saying 
that only the component of a normal to the level lines of the solution is relevant (Figure 5). 
In any analysis relating to events within the triangle T, a can be replaced by a T . There is 
an intuitive attraction to making a T as small as possible, that is, taking it in the direction 
Vu T . Then we have 

a ■ Vu T \ T fVu T 

Vu T • Vu T J U ~ Vu T • Vu T ' 
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A practical disadvantage to (9) is that in any region where the solution is almost uniform, 
the numerator and denominator are both very small, and numerical expressions based on 
(9) become ill-formed. This leads to slow, or halted, convergence. In [2], the expedient was 
taken of abandoning (9) gradually in favor of a in regions where jV^j is small. Below, we 
will argue that (9) is anyway not optimal; the improved version leads also to better formed 
expressions. 

We commence the analysis by noting that if (10) is substituted into (3) we have 

K — & T • Ei 


= a ■ + AaJ • n^. 

In the second inner product, aj is orthogonal to Vu and is orthogonal to Sj where 
is the vector along the side opposite Uj). Therefore 


K = a ■ Ui + \u ■ & 


so that 


h = K 4- A (u 3 — u 2 ) 


kj = K + ~ u 3) } 


( 12 ) 


(13) 


&3 = ^3 + A(u 3 — til) J 

Here, k*,^, are the original coefficients computed from (3), and equations (13) define 
a one-parameter family of equally valid alternatives. 

We pose the problem of choosing the k’s so as to maximize the allowable time-step of 
the scheme. This is clearly a desirable objective in itself but heuristically it seems to lead to 
gains in accuracy also. (In one dimension, this criterion would lead us to classical upwinding 
and in Section 3, we saw that it led to the optimum linear schemes in two dimensions.) From 
(7), the scheme is positive for Type II triangles, if 


At < — min 
n 



where i,j are the downwind vertices, and where n is the maximum number of triangles 
meeting at a point. Therefore, the relevant criterion is that 


K = max 


ki kj 

S’ 5 


(14) 


should be positive but as small as possible. 

Consider a diagram of which the axes are (&,/£), (kj/Sj). As A varies, a straight line 
locus will be swept out whose slope is 

_ Sjjuj - Uk) 

Sj{uj-u k y 


4 


Figure 6 shows such loci for given k*,kj (i.e., given geometry) for a range of A (i.e., a 
range of data). For each locus, a star marks the point that minimizes K . 

If the star falls on the vertical axis, the decision is to update only the point j ; on the 
horizontal axis we update only the point i. On the 45° line we have 


K 

Si 


S, 


so that, from (7) 


SuJ U{ — U k 

SuJ Uj - u k 


(15) 


(where SuJ means the change at vertex^ due to this triangle). This corresponds to keeping 
the direction of Vu constant during the update process. Figure 6 demonstrates that this 
is the best policy whenever A is negative, i.e., (u» — u k ), ( Uj — u k ) have the same sign. An 
equivalent criterion is that the level line of the solution through v k does not pass through 
the triangle. If the level line does pass through the triangle, we will only update one vertex. 
Which vertex is chosen will depend whether the locus in Figure 6 passes above or below the 
origin, and it is easy to demonstrate that this depends on the sign of <jy^ 

The update algorithm is in fact, in pseudo- Fortran 


If (only ki.gt. 0) 
update Vi 
goto end 
endif 

If (only k k .lt. 0) then 

If (( Ui — u k )(uj — u^.gt.O) then 
update Ui y Uj 

else if ((itj — u k )<f> T .lt. 0) then 
update Ui 
else 

update uj 
endif 
endif . 

There are actually twelve different paths through this logic, although only six different 
outcomes. It can be verified that the updates at the vertices depend C° continuously on all 
aspects of the data. The only part of the algorithm that is badly formed arises if both u, , u 
are to be updated, in accordance with (15). The explicit formulae are then 

SuJ = — r<f> T At 

Si(ui - u k ) + Sj(uj - u k y 

£ T — Uk . 

3 Si(uj - u k ) + Sj(uj - u k ) T 



(16a) 
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and it is possible for the denominator to be small. However, we only reach this situation if 
(iij — u*. )('u J — Ufc) are of the same sign, so they must be small individually, and (j>x must also 
be small. Thus, that triangle is in equilibrium, and can be ignored. In practice, we may test 
to see if <f>r is less than some moderate multiple of machine zero. 

We call this the NN scheme, signifying non-linear and narrow. 

5. Results and Discussion 

The schemes will be demonstrated on two test problems shown in Figure 7. 

The first problem is to solve the equation 

U t +U x + ~ Uy = 0 

on the unit square with u = 1 specified on the left edge and u = 0 on the bottom edge. This 
is the direction for which the IV, and NN, schemes will show least advantage. The second 
problem is to solve the equation 

u t + x u y — yu x = 0 

on the domain ne[— 1,1], ye[0, 1]. The characteristic lines here are circles centered on the 
origin. 

On all inflow boundaries we specify u = 0, except that on y = 0, we set u = 1 for 
xe [- §, — |]. The steady solution of this problem is therefore 

u = 1 if 1 < 9(x 2 + y 2 ) < 4 
: 0 otherwise. 

In particular, when we plot u(x, 0) we expect, or at least hope, to see the input data for 
x < 0 mirrored in the solution for x > 0. An interesting feature of examining the solution in 
that way is that the dissipation error has been effectively integrated over all flow directions. 

Figure 8(a) contains results from the regular upwind scheme for Problem 1 on a mesh of 
31 x 31 nodes. The solution is very diffused, but would be even worse at a direction of 45°. 
Figure 8(b) is obtained from the SI scheme, but the mesh is triangulated in no preferential 
direction (Figure 7(c)). The results are very slightly improved, although this cannot be 
determined from the plots. 

A very noticeable improvement comes from deploying the SI scheme on a mesh trian- 
gulated along the upgoing diagonals. The results are shown in Figure 8(c), and precisely 
the same results are obtained using the N scheme on a square mesh. The discontinuity is 
resolved over about half the number of mesh points. It is worth remarking that this reduces 
by a factor of about 16 the computing cost of obtaining a solution with given resolution, 
because half as many points are needed in each direction, half as much time is needed to 
expel transients, and twice the timestep can be taken. 

In Figure 8(d) we see results from the nonlinear scheme NN. Even on the isotropic mesh, 
these are the best so far, but by choosing the correct triangulation, the further improvement 
of Figure 8(e) is reached. Since this improves the resolution by a factor of 4, the cost at 


6 



given resolution is reduced by about a factor of 64. Recall that the results from the optimum 
scheme are actually worse for this flow angle. It is gratifying that the convergence history 
is scarcely affected by any of these improvements. The residual stays almost constant for 
about 80 timesteps, by which time the initial guess u = 0 has been purged by the wave 
crossing the domain; after a further 70 to 80 timesteps, the solution is effectively converged. 
All tests were made at a Courant number of 0.5, but in fact the optimum scheme could have 
been run faster. 

In Figure 9(a), the results are shown of the rotating advection problem on a mesh of 61 
cells by 31 cells. Only the results for an optimized mesh are shown; the squares are cut by 
their upgoing diagonals for x < 0, by their downgoing diagonals for x > 0. The square pulse 
is preserved quite remarkably for a first-order scheme. Indeed, one may suspect that the step 
function is some kind of eigenmode for the scheme, and that all data will be squared-off, as 
happens for certain one- dimensional limiter schemes. To check this point, a narrow Gaussian 
pulse was substituted as data and was again preserved rather well (Figure 9(b)). 

To make some kind of comparison with regular upwinding, the calculation in Figure 
9(c) was performed. To give the regular scheme a chance, four times as many mesh points 
were used in each direction, i.e., 241 x 121. Even so, a very inferior result is achieved. In 
comparison with the straight advection test, the optimum schemes show in a better light and 
the regular scheme in a worse light. Their resolving powers differ here by a factor of about 
10, and the cost of achieving given resolution changes by a factor of about one thousand. 

6. Concluding Remarks 

The schemes derived here have features with no close one- dimensional counterparts. They 
actively select the stencil to be used, and the connectivity of the mesh, and depend non- 
linearly on the data, all in the search for an optimal first-order scheme. They are of course 
to be viewed as the first phase in a quest for higher-order schemes that may exploit the same 
features. Nevertheless, their practical use in simulations of mixing and turbulent flows may 
not be out of the question. The resolution that they achieve compares favorably with that 
of formally higher-order schemes, and because all their operations are highly localized, the 
schemes are very suitable for parallel implementation. 
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Figure 7. The Test Problems. 
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